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Abstract 



As part of an ongoing effort to characterize the high temperature phase 
of QCD, in a numerical simulation using the staggered fermion scheme, we 
measure the quark baryon density in the vicinity of a fixed test quark at high 
temperature and compare it with similar measurements at low temperature 
and at the crossover temperature. We find an extremely weak correlation at 
high temperature, suggesting that small color singlet clusters are unimportant 
in the thermal ensemble. We also find that at T = 0.75 T c the total induced 
quark number shows a surprisingly large component attributable to baryonic 
screening. A companion simulation of a simple flux tube model produces 
similar results and also suggests a plausible phenomenological scenario: As the 
crossover temperature is approached from below, baryonic states proliferate. 
Above the crossover temperature the mean size of color singlet clusters grows 
explosively, resulting in an effective electrostatic deconfinement. 
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Typeset using REVT^X 



I. INTRODUCTION 



Numerical simulations of the quark plasma have suggested seemingly contradictory mod- 
els. While bulk thermodynamic quantities, such as the energy density [1] and baryon suscep- 
tibility yield values consistent with a nearly free gas of quarks and gluons, measurements 
of screening propagators, particularly, measurements of the wave functions of exchanged 
objects, are consistent with the confinement of color singlets B. Indeed, simulations and 
analytic work in the pure glue sector have demonstrated that space-like Wilson loops obey 
an area law in the high temperature phase, a signature of confinement M. 

One resolution of this seeming paradox describes the quark plasma as an ensemble of 
color singlet clusters of various sizes. Bulk thermodynamic quantities, such as the energy 
density, would receive contributions from all clusters, whereas long-range screening would 
be controlled by the lightest clusters. How large is the typical color singlet cluster? What is 
the typical spatial extent and quark and antiquark content? To answer these questions, it 
is necessary to seek observables that have not hitherto been studied in this context. Thus, 
we measured the distribution of induced quark charge (baryon number) in the vicinity of 
a fixed test quark, at low and high temperature, and at the crossover temperature. This 
observable has also been studied by the Vienna group in an effort to discern changes in the 
QCD vacuum induced by color charges ]5|]. At low temperature we expect that, as a result 
of confinement, a dynamical antiquark or, less often, a pair of quarks, screens the test charge 
at short distance. Thus, the induced dynamical quark number density should be large and 
negative close to the test charge. If screening is entirely due to a single antiquark, we should 
observe that the total induced quark number Q is —1. By contrast, if color singlet clusters 
are large either in size or in the number of quarks and antiquarks, we would expect only a 
small induced charge density near the source. 

In Sec. H we describe the observable in detail and in Sec. [HI] we present the results of 
the numerical simulation. Among the more striking results is the surprising weakness of the 
induced quark number density at high temperature. We also find that at temperatures near, 
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but below the crossover, the total induced charge is significantly different from what would 
be expected if only a single antiquark were responsible for screening. || 

To help understand these results we turned to a simple flux tube model of Patel 0. This 
model incorporates some of the essential features of QCD, but has the added appeal that 
it can be formulated simply, either in a field-diagonal basis analogous to the Wilson QCD 
action or in a phenomenologically suggestive flux-diagonal basis. The latter representation 
permits the study of the growth and complexity of the color singlet clusters of the model. 
The induced quark density and charge can also be studied in this model. Results of this study 
are presented in Sec. |IV|. We find results surprisingly similar to those for QCD. A direct 
examination of color singlet cluster size in this model suggests an explanation for the QCD 
results, including a description of the nature of the phase transition, and of the structure 
of the high temperature phase. In particular, it suggests that heating of the confined phase 
results in the appearance of quark clusters, including a surprisingly high number of baryons 
and antibaryons. As the crossover temperature is passed, these clusters grow explosively, 
both in size and in quark content, resulting in a suppression of baryonic correlations, a rise 
in the baryon susceptibility, and an effective electrostatic deconfinement. 

A concluding discussion is given in Sec. |V|. 

II. INDUCED BARYON DENSITY FOR STAGGERED FERMIONS 

Here we derive the staggered fermion observables to be measured in the numerical sim- 
ulation. 

A. Baryon density from the local chemical potential 

The construction of the local quark number density starts with the introduction of a 
baryon chemical potential in the standard way ||, but with a spatial dependence. Such a 
definition assures that the total baryon charge so defined is exactly conserved on the lattice. 
We start with the lattice action for staggered fermions in the notation of Ref. 0, modified 
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through the introduction of a local chemical potential. Let r = (r, t) and r' = (r', t') denote 
lattice coordinates. The action is 

S{U, i>, j,) = S g (U) + £ VS(r)M(U) rjr ^(r') (1) 

r.r' 

where S g (U) is the pure gauge action. The fermion flavor is implicitly summed over. The 
fermion matrix for a single flavor with a local chemical potential /x(r) is given by 

M(U) r y = 2ma5 r y + ^ \U r>v 8 r y- v — Ur-v,v^r,r'+u) + M t (U) r y (2) 

where the link matrices include the usual Dirac phase factor 

Ur,v = T]r,vUr,v (3) 

The Dirac phase factors r] r j also include the sign for the antiperiodic boundary condition. 
We consider two alternative formulations "static" and "slice" for the the time-hopping part 
of the fermion matrix M t . First 

Notice that the fugacity factor has been spread uniformly in the time dimension. We also 
consider an alternative "single-slice" definition that introduces the fugacity factor on a single 
time slice: 

MLe(U) r y = [U r , t e N ^5 ry _ { S t , fi - Ul_ tf - N ^5 ry+i 5 tfl ] (5) 
+ [U r . t 5 ry _ i (l - 6t>, ) - Ul_ t / ry+i {l - Stfl)] 

Of course, which formulation we choose depends on what we want to measure. If the chemical 
potential /x(r) is introduced on a single time slice as in Eq. (5), it is a source for the local 
baryon density at a single time. If it is instead spread over all time as in Eq. (|]), it is a 
source for the time-averaged (static) baryon density. When /i(r) is independent of r it can 
be shown that the determinant of the fermion matrix is the same in either case. Thus, for 
example, the baryon susceptibility is obtained by differentiating the free energy twice 
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with respect to such a constant chemical potential, so it is the same in either formulation. 
Moreover, for any local time-independent observable, such as the baryon density in the 
presence of a static charge, the expectation value is the same. A difference appears in the 
dynamical baryon density-density correlation, which may depend on the time separation of 
the operators. The single-slice formulation gives the equal-time density-density correlation, 
and the static formulation gives the correlation of the time-averaged densities. For those 
observables that are independent of the formulation, we find that the static version has the 
practical advantage that, because of time averaging, it produces a higher signal to noise 
ratio. 

The baryon density per unit lattice cell a 3 at zero chemical potential at a chosen spatial 
coordinate r is given in terms of the individual flavor densities by 

N f 

P(r) = EP*M (6) 

i=i 

where 

Pl (r)=p- 1 d\ogZ/d^(r)\ Mr)=0 . (7) 
Here f3 = l/T. Now, for two flavors of staggered fermions the partition function is 

Z = J [dU] exp[-S g (y)][dest M U (U, pi u ) det M d (U, /i d )] 1/4 (8) 

and we have 

( PM (r)> = (4/T 1 (Tr[M- 1 d dM u , d /d f i u , d (r)]) u \, u d{r)=0 (9) 

The result is 

W')) = (V4)0. (itclM^),^^);,] + ItclM^.^^.^^ (10) 

where r is summed over and Tr c denotes a trace over color indices only. As a device for 
treating both cases "static" and "slice" for the time-hopping term in the same expression, 
we have introduced the weight 
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0static,r = 1/iVt 0slice,r = 5 T ,Q- (H) 

Because the hopping term in the staggered fermion matrix has the time-reversal symmetry 

K,r> = (-Y~ r 'M r ~!r (12) 

the two terms on the rhs of the expression for the density are negative complex conjugates 
of each other. We get 

( Pu , d (r)) = (z/2)Im(P dyn (r)) (13) 

where 



-Pdyn(r) = Tr c 



Clearly for the particular time-independent case of Eq. ( |i"3"D both the "single slice" and 
"static" formulations yield the same result. In particular, since the expectation value on the 
rhs is related to the Polyakov loop expectation value, which is real for the usual action at 
zero chemical potential, we get zero for the baryon density at zero chemical potential, as we 
should. 

B. Density in the presence of a test quark 

Now we want to consider the correlation between a point test quark at the origin (say) 
and the baryon density at r. Introducing the test quark simply involves modifying the action 
by including a Polyakov loop factor Pfi xe d(0), defined through 



-Pfixed( r ) — Tr c 



N t 



II u (r,ty,t 

t=0 



(15) 



We then have the partition function for the ensemble with the test quark: 

Z q = J [rff/]exp[-S , ,(f/)][detM u (f/, /Utt )detM d (f/, /Ud )] 1 / 4 P fixcd (0) (16) 

The baryon density p q (r) in the presence of the fixed quark can be calculated on the original 
ensemble Z through 



p,(r) = (iN f /2) (P fixcd (0)ImP dyn (r)) [/ / (P &xcd (0))u ■ ( 17 ) 

Now the expectation values are taken with respect to the original test-charge-free ensemble. 
The correlation between RePfi xed and ImP dyn vanishes because of the complex conjugation 
symmetry of the integral over the gauge variables. What survives is the correlation of the 
imaginary parts divided by the expectation of the real part of the Polyakov loop, namely 

p q (r) = -(N f /2) (ImP fixcd (0)ImP dyn (r)) f/ / (ReP fixed (0)) [/ . (18) 

C. Density-density correlation 

The connected density-density correlation for flavors i and j is given by 

^r) = {ft(r) ft (0))-(ft(0)>^(0)) 

= /T 2 c> 2 log Z/d/^r, 0)0^(0, 0)U (r)=0 (19) 

In terms of the integrated correlation 

ptot,ij = J d 3 rpij(r) (20) 

the baryon singlet and nonsinglet susceptibilities [§] are 

XS = P(Ptot,uu + Ptot,dd + Ptot,ud + Ptot,du) (21) 
XNS = P(Ptot,uu + Ptot,dd — Ptot,ud — Ptot,du) (22) 

In terms of the fermion matrices for the separate flavors the correlation receives four contri- 
butions 

Pij( T ) — Sllij + S2Uj + S%2ij — ^disc,^'; (23) 

The indices nm in S nm ij count the number of explicit coordinate points n and number of 
color traces m. 
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S niJ = (2(3)-% (Tr[M^d 2 M l /d^(r)d^(0)]) u L, d M=o (24) 

Sna = (Tr[Mr 1 dM i /d f x i (r)Mr l dM i /df M (0)}) u \, uAr)=0 (25) 

Sv* = (4/3)- 2 (Tr[Mr l dM i /d f , i (r)]Tr[M- 1 dM j /d f , j (0)}) u \, uAr)=0 (26) 

S discij = <pi(0)> (^(0)) . (27) 

Carrying out the derivatives and using the time-reversal property fljjD gives 

S 1Uj = 2~ 1 6 lJ 6 r , (ReP dyn (0)> (28) 

-2"%^ (ReTr JM ( -* T+ (29) 

S 2 2ij = -4" 1 (ImP dyn (r)ImP dyn (0)) l/ (30) 



The two-point single-trace term S2iij is a hadron propagator with a source at (0, r) and a 
sink at (r, r'). The source and sink are both just the point-split baryon density operator. 
To evaluate this term we require the quark propagator from the source U(o, T y,t at (0, r) and 
the antiquark propagator from a point source at (0, r + 1). These propagators are combined 
in two ways at (r, r') to complete the evaluation of the two contributions. 

To evaluate the two-point, two-trace term, S 2 2ij we use the random source trick of Ref. [0 . 
We introduce a set of n rSLnd independent complex Gaussian random SU(3) vectors P^(r, r) 
I = 1, . . . , n ra nd on the time slice r = 1 for the "slice" form of the action and for all r for the 
"static" form. Then on a given gauge configuration U we generate the Fourier transform of 
the estimate of ImP dyn (r). 

J dyn/ (k) = £exp(ir • k)0 T Im[P^(r)*M ( - 1 T+1)j(riT) C/ (riT)l tP,(r)]. (31) 

r 

(Note that the imaginary part is taken before carrying out the Fourier transform.) 
Finally, we estimate the Fourier transform of S 2 2ij from 

Wk) = - T- ,] —r, E ( W k ) hyn^)) v ■ (32) 

^ ' <rand \ ' ^rand *-) g-^gi 

The single-point, single-trace term can be generated trivially from the random source method 
through the estimate 
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Smj = E E T)*M ( -V fl)i(r>T) ^ lT) , t i?£(r > r)]. (33) 

^ "rand ' / r 

Thus the computation of Smj and S^ij starts with the evaluation of the quark propagator 
from the same parallel-transported random source, namely, tJ/ rtT ^ t Rk(T, r). The same Fourier 
transform Idyn,i is used to obtain the correlation with the static source. However, evidently 
the computation of the hadron propagator S 2 uj must be done with point sources — not 
random sources. 



III. RESULTS OF SIMULATIONS 

Simulations were carried out at fixed (3 = 5.445 and quark mass am q = 0.025 for two 
flavors of staggered fermions on lattices of size 16 3 x N t , where N t = 8,6,4. This choice 
of lattice parameters corresponds to the crossover temperature at Nt — 6 @. Thus, the 
simulations are done at three temperatures T = 0.75 T c , T « T c , and T = 1.5 T c , respectively, 
at the same lattice scale, making it meaningful to superimpose plots of baryon density vs 
distance from the source. Spectroscopic simulations at the same temperature [TO] allow us 



to set the scale, viz. T c = 145 MeV and a = 0.227 fm. Simulations were also carried out at 
Nt — 4 with (3 = 5.15, 5.22, 5.25 and 5.29 (the crossover) to provide an independent check 
of trends in the total induced quark number. Table | shows the extent of the simulation 
sample. 

Figure [l| summarizes our results for the induced quark number density at these three 
temperatures. Particularly striking is the dramatic decrease in the correlation at high tem- 
perature. Thus, we see no evidence for small color singlet clusters in the high temperature 
plasma. 

The total induced quark number normalized to one for a single quark was computed in 



two ways: first by a direct integration of the density (18) and second by fitting the density 
distribution to the functional form 
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(the sum of free lattice propagators for scalar fields of mass rrij) from which the total induced 
quark number is 

Q = N^^/m 2 . (35) 

3 

At most two mass terms were used. All fits started at zero radius. The resulting total quark 
number determined from the two methods was in each case consistent within errors. The 
value quoted is the one with the smaller standard deviation. It was found that except for the 
high temperature points (6/g 2 = 5.445, N t = 4 and 6/g 2 = 5.29, N t = 6) direct summation 
over the entire volume gave poorer statistics than fitting, because direct summation suffers 
from statistical noise introduced by contributions far from the fixed charge that should sum 
to zero. By fitting to a functional form that falls to zero at infinity, we controlled this noise. 

The resulting fitted curves are plotted in Fig. |l|. The total quark number values are also 
given in the legend and in Table 0. At the high temperature point the total induced quark 
number is nearly two orders of magnitude smaller that the quark number at low temperature. 
At low temperature we expect that the test charge is attached to a color singlet cluster. A 
single antiquark would contribute —1 to the total induced quark number, and a pair of 
quarks forming a baryon, +2. A thermodynamic mixture of these two configurations would 
give an intermediate value. Based on the error ellipse for a two-parameter fit to the lowest- 
temperature density, at the two standard deviation level, we find that the total induced 
quark number is greater than —0.55, significantly different from —1. 

As a check of this result, we also carried out a series of simulations at varying (3 with 
Nt = 4. At such a strong coupling the Polyakov loop expectation value is large, leading to 
a stronger signal in the correlation. Results for the induced quark number are plotted in 
Fig. |2] and listed in Table [TIT]. The low temperature values are somewhat larger in magnitude 
that in the 6/g 2 = 5.445 simulation, but still show a significant departure from —1. 

For the sake of comparison let us estimate the contribution to the induced charge from 
the lowest S-wave mesonic and baryonic screening clusters in the ensemble. These clusters 
are obtained by replacing one quark in the n, p, N, and A by a fixed spinless, flavorless color 
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triplet quark. The result is a modified J = 1/2, 1 = 1/2 meson with a four-fold multiplicity 
g p i = 4, a modified J = 0, / = "nucleon" N' and a modified J = 1, I = 1 delta A'. In the 
continuum limit the degeneracies are g^i = 1 and #a' = 9. However, at finite lattice spacing 
the A' is split, owing to the discrete nature of the internal symmetries in the staggered 
fermion scheme. We have not measured masses of the modified states, but have masses for 



their light quark counterparts [T(| for the same coupling 6/g 2 = 5.445, namely, am p = 0.918, 
om,, = 0.4488, am^ = 1.375, and owa = 1.43. If we assume that the splitting of the p 
and 7r is entirely due to the color hyperfine interaction, then we estimate the mass of the 
modified meson to be m p i = M + (3m p + m 7r )/4 = M + 0.80/a, where M represents the 
contribution from the point charge. Similarly, we have m^i = M + = M + 1.375/a and 
rriA' = M + (2mA + J7ijv)/3 = M + 1.41/ a. To be conservative, let us assume that the A' is 
fully degenerate. As noted before, the induced quark number is 

Q = -Prn + tyb, (36) 

where p m and pi, are the probabilities of screening via the mesonic and baryonic clusters. 
For this estimate we take p m +Pb = 1- The probabilities are estimated from the Boltzmann 
weights: 

g A , e - m ±'/ T + g N ,e- m x'' T _ 9e~ 1A1 / aT + e ~ l -™' aT 

Vm ' Pb ~ g P >e~ m P ,/T ~ 4 e -0.80/aT ( 37 ) 

The unknown regularization-dependent energy M has cancelled in the ratio. The resulting 
estimates are Q = —0.94 at T = 0.75 T c = l/8a and Q = —0.81 at T = T c = 1/Qa. These 
values are considerably lower than were found in the simulation. To bring the estimates 
into closer agreement would require adding more baryonic states. Thus the full simulation 
suggests that already at a temperature of 0.75 T c , there is significant baryonic screening of 
the fixed charge. 

As we have remarked fl22"|) , the integral of the self-correlation of the dynamical quark 
number density gives the baryon susceptibility. Results for the susceptibility on 8 3 x 4 
lattices with ma = 0.025 were reported in Ref. f2j. The high temperature values found are 
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consistent with what would be expected for an ideal quark gas, given the very large effects 
of the nonzero lattice spacing. As a check, we compare our new results on 16 3 x 4. Figures 
|3] and ^ show the comparison for the N t = 4 series. They are apparently consistent. Results 
for the N t = 4 series are also given in Table [TTT] and for the Q/g 2 = 5.445 series, in Table [TI|. 
Not surprisingly, the static form of the density operator gives a better signal than the single 
time slice form, because it involves an average over all time links. Values are not available 
for (3 = 5.29, which was run before improvements in the code incorporated the static form 
and gave an acceptable signal to noise ratio for the slice form. 

The correlation between a test charge and the scalar density (ijrtp^ has been measured 
by the Vienna group ||. Our results, shown in Fig. [|, are consistent with theirs. 



IV. FLUX TUBE MODEL 
A. The model 

Some years ago Patel J7J proposed a flux tube version of the three-state three-dimensional 
Potts model to explain the mechanism of the deconfining phase transition in QCD. In this 
model, each site r of a cubic lattice holds either a quark, antiquark, or none at all, and each 
link £ r)fl , a triplet or antitriplet flux, or none at all. That is the quark number n r and the 
flux £ r )At take on values {—1, 0, 1}. Flux is conserved modulo 3. 

3 

+ 4 - n T = mod 3 (38) 

where £ r ,-n — —^r-p,,n- The hamiltonian is given in terms of the quark mass m and the 
string link energy a by 

# = 5>|4, M |+X>H- (39) 

r,fj, r 

The partition function is then 

Z(/3)= £ exp(-PH) (40) 

where the prime signifies a sum constrained by (j38|). 
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B. Equivalence to the Potts model 



This model is equivalent to the three-state three-dimensional Potts model. The equiv 
alence can be seen by replacing the Gauss' law constraint by 

1 



^o = ^ E z 

6 z£Z(3) 



(41) 



in modulo three arithmetic. Here Z(3) = {1, e ±27r / 3 }. Introducing this identity on each site 
with the summation variable z r allows us to rewrite the partition function as 



Z{p)= E exp -/? 



EH4J + E 



m\n r 



r,/i 



n 



(42) 



The unconstrained sums over links and quark numbers can be carried out as follows: 

Eexp(-/3|4J)(^ r * +/i )^ = l + 2Re(z r < +A )exp(-pV) 
Eexp(-/3m|n r |);z~ nr = 1 + 2Rez r . 



(43) 
(44) 



Thus the partition function reduces to a product of polynomials in the Z(3) variables z r . A 
hamiltonian can be constructed through the identity over Z(3) 



log(l + cz + cz*) = a + bz + bz* 



(45) 



where 



exp(2a) = (l + 2c)(l - c) 2 
exp(36) = (1 + 2c)/(l - c) 



(46) 



Thus if we define 



3 y 1 — exp(— per) 

h0 = 2 -J l + 2e ^^ 



3 y 1 — exp(— (3m) 



(47) 
(48) 



then we have, up to a constant 
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where 



Z{(3) = £ exp(-/3'#') 

{*■} 



(49) 



H' = - JRe ( Z rZ*+fi) - £ hReZ *- ( 50 ) 



This expression is recognized as the hamiltonian of the three-state Potts model in three 
dimensions with a coupling J and a magnetic field h coupled to the real part of the spin. 

Thus a low quark mass corresponds to a high Potts magnetic field and a low flux-model 
temperature corresponds to a high Potts temperature. At zero field h a first order phase 
transition is found in this model at J/3' ~ 0.367. The phase transition persists for a small 
magnetic field h(3' < 0.002, but is not evident in numerical simulations for larger values 
of the field ||11||. These parameter values can be converted to the flux tube model values 



through the inverse of Eq. (|48|) : 

* = ln Udw-ij (51) 

Pm^J^f^ 2 ), (52, 

giving a phase transition along a curve starting at about [3a = 1.63, (3m = oo to about 
(3m > 6.9 or mja > 4.2. 



C. Relationship to QCD 

Patel proposed using this model as a paradigm for the QCD phase transition. The Potts 
model in its more conventional form was also offered some years ago as a model of the 
deconfining phase transition [|12j. The latter formulation is obtained from a high-quark- 
mass, strong-coupling, high-temperature, Z(3)-restricted approximation to the conventional 
field- diagonal Wilson action, with the Potts spin corresponding to the Polyakov loop. On 
the other hand the flux-tube formulation of the Potts model corresponds to an alternate 
representation of the Wilson action in the charge-and-flux-diagonal basis. Because of the 
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combinatoric complexities of linking SU(3) charges and fluxes to form color singlets, such a 
basis for the Wilson action never received wide attention. However, in the simple Z(3) basis 
of the flux tube model the combinatorics become trivial. Moreover, in the flux-tube form 
the model offers the highly suggestive possibility of studying the size and structure of color 
singlet clusters as a function of temperature. Its chief drawback is that, because it treats 
quarks as static objects, it does not incorporate chiral symmetry. 

At low temperature only small color singlet clusters populate the Gibbs ensemble. As 
the temperature is increased, clusters of increasing size occur. Eventually clusters connect 
to fill nearly the entire spatial volume. For heavy quark masses this phenomenon leads to 
a first-order deconfinement phase transition. For light quarks, cluster growth is somewhat 
inhibited, since pair formation breaks the flux links. It is found that there is no phase transi- 
tion. One is tempted to think of a percolation phase transition mediated by the connectivity 
of the flux tubes, but there is nothing to percolate: there is no current in the model to flow 
between linked sites and establish long-range order. The first-order phase transition occur- 
ring only at large quark mass (low magnetic field) resembles more appropriately a liquid-gas 
phase transition. 

D. Simulation and Observables 

To simulate the model in the flux tube basis, we used a Metropolis algorithm, with moves 
designed to preserve Gauss' law. We considered two types of elementary moves: the addition 
(or removal) of the lightest meson (flux link with quark and antiquark at the ends) and the 
addition (or removal) of the lightest "glueball" (four flux links directed around a plaquette). 
Adding was done in the literal sense: adding a quark to a site means increasing the quark 
number by one unit, modulo three, etc. Thus, the actual impact of these elementary moves 
depends on the configuration. The moves could result in shortening, lengthening, breaking, 
joining, or displacing a series of flux links. Notice that two or three such mesons can form a 
baryon, if they are added with the antiquark on the same site but with links and quarks on 
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unique sites, so we didn't create baryons through a separate move. Although by the same 
token four such mesons can form a glueball, since meson formation is suppressed at high 
quark mass, we kept both moves to allow a more efficient mass-independent coverage of the 
ensemble. 

A fixed charge is introduced at the origin by starting from a modified vacuum configu- 
ration in which the dynamical quark charge at the origin is set to —1 and all other charges 
and fluxes are initially zero. Observables include these: 

| n | mean number of quarks plus antiquarks per site 

\i\ mean number of links (either sign) per link 

n vtx mean number of three-point flux vertices 

p(r) induced quark number density 

Q the total induced charge (quarks minus antiquarks) 

(B 2 ) mean square baryon number (including test charge) 

iV the size of the cluster connected to the origin 

E. Results 

Measurements were made on a 10 3 lattice for a variety of (3 at m = 1.0 (10,000 sweeps 
for each (3 value) and m = 3.5 (1,000 sweeps for each (3 value). We use units in which a — 1. 
As we have noted above, both quark masses are in a region where a phase transition does 
not occur. The higher mass series comes closer to the critical point. With such a small 
volume we would notice a significant finite-size rounding of the first order phase transition 
that occurs at higher quark masses. However, we are primarily interested in the qualitative 
behavior of the model at smaller quark mass, since it corresponds more closely to the light 
quark QCD simulation. Thus the small volume suffices. 

A popular indicator of the phase transition in QCD is the Polyakov loop ReP. The Potts 
model analog is the magnetization, 

(Re, - (53, 
17 



From the map Eq ( fHf) onto flux tube variables, we see that the flux tube analog of this 
order parameter is the mean quark count (quarks plus antiquarks) 

/i i\ d\ogZ 

(W) = Wd^ (54) 

Figures ^| and [7] show the quark number \n\ vs ft for the two quark masses. The higher 
mass series comes closer to the critical point, leading to a sharper crossover. The crossover 
locations are determined from the inflection point of the curves (peak in the corresponding 
susceptibility) to be 1/T C — /3 C — 1.88(1) for m = 1.0 and (3 C = 1.662(2) for m = 3.5. These 
values are used to determine the ratio T/T c in Tables |V| and [V|. 

Shown in Fig. || and Tables |TV| and [V] is the total induced charge Q as a function of 



/3 in the presence of a test charge for the two quark masses. The results bear a striking 
resemblance to those of the QCD simulation. At high temperature, the induced charge is 
very small for both quark masses. For the heavier quark mass, the total induced charge 
Q = —0.75(7) at the highest beta (T = 0.92 T c ) shows a significant departure from —1. For 
the lighter quark mass at the highest f3 (T = 0.78 T c ) we also see a significant departure 
with Q = -0.69(2). 

Thus we find evidence of important baryonic screening in this model at temperatures 
close to the crossover. Is this surprising? Consider the relative Boltzmann weights for the 
lightest baryon-like and meson-like clusters attached to the test quark. They are shown in 
Fig. |9|. The masses are m m = M + m + a and m;, = M + 2m + 2a for the meson-like 
and baryon-like cluster, respectively. Here M stands for the mass of the fixed charge. For 
m = 1 and a = 1, corresponding to our the lighter quark mass, these are m m = M + 2 and 
rrib = M + 4. The corresponding multiplicities are g m = 6 and gj, = 30. At (3 = 2.4 we have, 
in the notation of Eq (|37|) 

Pb/Pm = 30e" 4/T /6e' 2/T = 0.04. (55) 

Thus, if only these two states played a role in low temperature screening we should find 
a total induced charge Q = —0.88, significantly different from what is observed. Clearly, 



still more baryonic states contribute. If we consider higher excitations, the number of dis- 
tinct baryonic states (i.e. states with total baryon number 5^0) grows very rapidly with 
mass — much faster than mesonic states. Thus the entropy of baryon excitation increases the 
importance of baryonic screening. This point can be dramatized by a direct examination 
of configurations. Shown in Fig. |ITJ is a representative configuration obtained at the mass 
and inverse temperature in question: m = 1 and (3 = 2.4. It contains nine qq mesons, one 
qqq baryon, and one qqq antibaryon, a highly improbable occurrence in the naive two-state 
model. 

Turning now to the induced density vs distance, we show results for the simulation at 
the lighter quark mass m = 1.0 in Fig. [11]. Again we see a resemblance to results of the 
QCD simulation shown in Fig. |l|. The total charge Q(r < 3) given in the legend is found 
by integrating the density to an arbitrary cutoff distance r < 3. The correlation is stronger 
at low temperature than in the QCD simulation, reflecting a shorter correlation length in 
lattice units, which could be adjusted by a change of scale. From a fit to a single-pole lattice 
Yukawa form, we find an effective mass of 2.9(7) in the flux tube model at (3 = 2.4 and 
m = 1.0, to be compared with 1.7(6) for the corresponding low temperature curve for SU(3) 
(Fig- 0). 

Further evidence of the importance of larger hadronic clusters can be obtained from a 
measurement of the mean size of the cluster attached to the origin. This size is defined as 
the number of sites connected through flux links to the origin. Our results are shown in 
Fig. [12] Notice that in this 10 3 lattice the mean cluster size grows dramatically, already 
filling 43% of the volume at T « 1.2 T c . Despite appearances, however, correlation lengths 
are nonetheless finite. 

Another indicator of clustering is the number of three-point flux vertices n vtx . This 
statistic is the sum of one third the total flux entering each site, if the total flux count at the 
site is positive. With coordination number six a site can contribute only 0, 1, or 2 to this 
statistic. As can be seen from Tables [TV] and |V| this number grows with temperature in much 
the same way for both the smaller and larger quark mass. On the other hand, a striking 
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difference is seen in the behavior of the baryon susceptibility (B 2 ) /V, as might be expected. 
For the lighter quark mass the fluctuation in baryon number is considerably larger at high 
temperature. The expected suppression of quark content with increasing quark mass is also 
evident in the quark number per site \n\. An expected consequence of decreasing quark mass 
is a decrease in cluster size, since pair creation breaks flux links. This effect is evident in 
the mean size of the cluster attached to the origin. At m — 3.5 and (3 = 1.4, corresponding 
to T = 1.19 T c , this mean size 650(3). When the quark mass is decreased to m — 1.0 at a 
comparable temperature (j3 = 1.6), the mean size decreases to 430(3). Thus, with decreasing 
quark mass the quark content of the clusters increases, and the size decreases somewhat. 



V. CONCLUSIONS 

In numerical simulations we have measured the quark number density in the vicinity of a 
test quark as a function of temperature. A strong correlation is found in the low temperature 
phase, but it is vastly reduced in the high temperature phase. We have also found evidence 
for a significant proliferation of baryonic clusters as the crossover temperature is approached 
from below. A companion simulation of the flux tube model has similar behavior and 
suggests an explanation. Heating past the crossover in this model results in an explosive 
growth of color singlet clusters. Thus at high temperature the addition of a single test 
quark has little effect on the ensemble, leading to an extremely weak correlation. We have 
an effective electrostatic deconfinement without a phase transition. 

We also find that in the flux tube model baryonic clusters proliferate as the temperature 
rises through T c , permitting more frequent baryonic screening of a test charge, suggesting an 
explanation for an apparent superabundance of baryons in full QCD at these temperatures. 
It is interesting to speculate that such a superabundance, particularly of antibaryons, should 
they survive final state interactions, provides an experimental signal for the crossover to the 



quark-plasma regime [13 



To be sure the flux tube model omits many features of QCD. It lacks dynamics, describing 
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only electrostatics. It also ignores chiral symmetry. Completely omitted are the important 
magnetic interactions that give rise to confinement in spacelike propagation. It would be 
useful to find an elaboration of the model more closely relevant to QCD. Nonetheless, it 
is highly suggestive both for further exploration of QCD and for the phenomenology of the 
quark plasma. 
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FIGURES 

FIG. 1. Quark number density induced by a fixed quark at the origin as a function of distance 
from the origin at three temperatures. Curves are fits to a single screening mass. The total induced 
quark number Q is also shown. 

FIG. 2. Total induced quark number vs 6/g 2 at Nt = 4. 

FIG. 3. Nonsinglet quark number susceptibility vs 6/g 2 at Nt = 4 with bare quark mass 
ma = 0.025, compared with Ref. || (crosses). Two forms of the quark number density are used: 
"static" (squares) and "slice" (circles). 

FIG. 4. Singlet quark number susceptibility vs Q/g 2 at Nt = 4 with bare quark mass 
ma = 0.025. Symbols are as in the previous figure. 

FIG. 5. Scalar density correlation \ij)ip) vs distance from a point charge, normalized to one at 
infinite distance. Error bars for the two lowest temperature points are as small or smaller than the 
plot symbol. For the sake of clarity, for the highest temperature only two typical error bars are 
shown. 

FIG. 6. Total number of quarks and antiquarks vs inverse temperature in the flux tube model 
for quark mass m = 1.0. The vertical line indicates the crossover. 

FIG. 7. Total number of quarks and antiquarks vs inverse temperature in the flux tube model 
for quark mass m = 3.5. 

FIG. 8. Induced charge vs temperature in the flux tube model. 

FIG. 9. Lightest meson-like and baryon-like clusters (two types) attached to a fixed quark 
(burst) in the flux tube model. 

FIG. 10. Typical flux tube lattice at T ~ 0.78 T c . The origin, indicated by the burst, has been 
displaced for clarity. 

FIG. 11. Flux tube model with m = 1.0. Induced quark number density vs distance from the 
origin. The legend Q(r < 3) gives the integral out to r = 3. 

FIG. 12. Mean size of cluster connected to the origin vs inverse temperature for the flux tube 
model. The errors are smaller than the plot symbols. 
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TABLES 



TABLE I. QCD simulation sample size (molecular dynamics time units) 



(3 N t time 

5.15 4 994.5 

5.22 4 856 

5.25 4 430 

5.29 4 500 

5.445 4 611 

5.445 6 1141 

5.445 8 2665 



TABLE II. QCD simulation results for 6/g 2 = 5.445 









static 


slice 


N t T/T c Q 


ReP 




Xs 


Xns 


Xs 


Xns 


4 1.5 -0.380(100) 0.6480(20) 


0.0768(5) 


0.233(10) 


0.2410(80) 


0.240(20) 


0.239(13) 


6 1.0 -0.130(60) 


0.1070(40) 


0.3150(30) 


0.030(7) 


0.0360(40) 


0.035(16) 


0.042(13) 


8 0.75 -0.006(3) 


0.0065(5) 


0.1860(5) 


0.009(5) 


0.0019(6) 


0.010(12) 


0.008(10) 




TABLE III. 


QCD simulation results for N t = 4 












static 


slice 


6/g 2 Q 


ReP 


<^> 


Xs 


Xns 


Xs 


Xns 


5.15 -0.790(70) 


0.0486(9) 


0.4945(10) 


0.000(20) 


0.025(6) 


0.05(3) 


0.040(13) 


5.22 -0.670(60) 


0.0680(20) 


0.4477(13) 


0.000(30) 


0.000(1) 


0.01(2) 


0.023(15) 


5.25 -0.600(60) 


0.0860(20) 


0.4110(20) 


0.070(20) 


0.045(9) 


-0.01(3) 


0.039(18) 


5.29 -0.046(6) 


0.4060(20) 


0.1990(30) 








0.171(14) 


5.445 -0.006(3) 


0.6480(20) 


0.0768(5) 


0.233(10) 


0.241(8) 


0.24(2) 


0.239(13) 
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TABLE IV. Flux tube model results for m = 1.0 



p 


T/T c 


n 




n vtx 


Q 


(B 2 ) 


N 


1.0 


1.88 


0.42220(20) 


0.42037(9) 


157.550(80) 


0.00(20) 


46.50(60) 


894(2) 


1.4 


1.34 


0.31693(15) 


0.30492(10) 


105.580(70) 


0.00(20) 


35.00(50) 


676(2) 


1.6 


1.17 


0.25783(15) 


0.23331(11) 


73.680(80) 


-0.01(15) 


28.40(40) 


430(3) 


1.8 


1.04 


0.18532(15) 


0.14651(12) 


39.080(60) 


0.03(13) 


19.50(30) 


73.5(1.0) 


2.0 


0.94 


0.10102(14) 


0.06091(10) 


11.980(40) 


-0.44(9) 


9.53(13) 


7.88(10) 


2.2 


0.85 


0.04315(10) 


0.01909(5) 


2.490(20) 


-0.51(5) 


2.75(4) 


3.25(3) 


2.4 


0.78 


0.01957(7) 


0.00714(3) 


0.636(8) 


-0.69(3) 


0.834(15) 


2.26(2) 


TABLE V. Flux tube model results for m = 3.5 


P 


T/T c 


n| 




n vtx 


Q 


(B 2 ) 


A^o 


1.40 


1.19 


0.0530(8) 


0.2687(4) 


71.1(3) 


0.10(20) 


1.77(6) 


650(3) 


1.50 


1.11 


0.0427(7) 


0.2149(5) 


51.3(2) 


0.00(11) 


1.35(7) 


535(3) 


1.60 


1.04 


0.0326(7) 


0.1386(15) 


27.9(4) 


0.20(10) 


1.09(7) 


337(4) 


1.62 


1.03 


0.0303(6) 


0.1202(15) 


22.7(4) 


0.00(20) 


0.97(6) 


289(5) 


1.64 


1.01 


0.0282(7) 


0.0950(20) 


16.5(5) 


0.20(20) 


0.94(6) 


206(7) 


1.66 


1.00 


0.0242(5) 


0.0670(20) 


10.1(4) 


-0.30(15) 


0.71(5) 


117(5) 


1.68 


0.99 


0.0183(7) 


0.0450(20) 


5.8(5) 


-0.17(9) 


0.54(5) 


59(6) 


1.70 


0.98 


0.0169(4) 


0.0341(11) 


3.8(2) 


-0.26(9) 


0.42(3) 


35(3) 


1.80 


0.92 


0.0099(3) 


0.0137(3) 


0.89(5) 


-0.75(7) 


0.12(3) 


8.2(4) 
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